%--------------------------------------------------------------------------
% JOINT DEBT & EQUITY PRICING
%--------------------------------------------------------------------------


% initialize
dumI            = zeros(nS,nK);
dumI(k'>k_iss)  = 1;
lamDx           = lamD - 1./k';
equity_iss      = k' < 1;

% eps_it : Idiosyncratic shock, corr with eps_ct is rho (under P and Q)
[Cheb_nodes,~]  = qnwcheb(nE,-1,1);
a               = -5;
b               = 5;
Eshocks         = (b-a)*(Cheb_nodes'+1)/2 + a; % Heer and Maussner p.600
temp            = pi*(b-a)/(2*nE);
weights         = temp.*sqrt(1-Cheb_nodes'.^2);
probE           = normpdf(Eshocks).*weights;
probE           = probE';
Eshocks         = Eshocks';

% Asset-to-earnings ratio
X       = (exp(mu_XQ + sig_X.*Eshocks')*probE) ./ Rf(:,1); % [nS,1]
lamA    = (1-tau_c)*(1-tau_d)*((eye(nS)-X.*Q)\ones(nS,1)); % [nS,1]

% next period's earnings growth and coverage ratio under Q
exp_gp_Q       = exp(mu_XQ(mgrid) + sig_X(mgrid)*Eshocks'); % [ nK*nS , nE ]
kp_Q           = kgrid.*exp_gp_Q; % [ nK*nS , nE ]
exp_gp_Q_iss   = exp(mu_XQ + sig_X*Eshocks'); % [ nS , nE ]

% initialize
Z_S     = NaN(nK*nS,nS);
Z_S_iss = NaN(nS,nS);
Z_D     = NaN(nK*nS,nS);
Z_D_iss = NaN(nS,nS);
Z_F     = NaN(nK*nS,nS);

iter    = 0;
max_iter = 500;

while iter < max_iter
    iter = iter + 1;
    
    kp_Q_iss  = k_bar.*exp_gp_Q_iss; % [ nS , nE ] today's state, tomorrow's E shock
    
    for loop = 1:5
        
        % iterate on equity
        for iter_S = 1:5
            lam_Dx_kbar = NaN(nS,1);
            for s=1:nS
                fun_lamS       = griddedInterpolant(k,lamS(s,:),meth,ex_meth); % here 's' is tomorrow's state
                lamSp          = fun_lamS(kp_Q);
                Z_S(:,s)       = (exp_gp_Q.*lamSp)*probE;
                lamSp          = fun_lamS(kp_Q_iss);
                Z_S_iss(:,s)   = (exp_gp_Q_iss.*lamSp)*probE;
                fun_lamDx      = griddedInterpolant(k,lamDx(s,:),meth,ex_meth); % here 's' is today's state
                lam_Dx_kbar(s) = fun_lamDx(k_bar(s));
            end
            EMS     = reshape( sum(Z_S .* Q_big,2) , [nS,nK] )./Rf1_mat; % [nS,nK]
            EMS_iss =          sum(Z_S_iss .* Q,2)            ./Rf(:,1); % [nS,1]
            lamSnew = max(0, (1-dumI).*( (1-tau_c)*(1-tau_d)*(1-1./k').*(1-equity_iss) + (1-tau_c+xi_e)*(1-1./k').*equity_iss + EMS ) ...
                + dumI .*( (1-tau_c)*(1-tau_d)*(1-1./k') + (1-xi_d-k_bar./k').*lam_Dx_kbar + EMS_iss ));
            dS      = max(abs(lamSnew(:)-lamS(:)));
            lamS    = lamSnew;
        end
        
        % dummy for state-dependent default threshold
        kd_s    = sum(lamS==0,2); % index of k-value that equals default threshold
        dumD    = zeros(nS,nK);
        for s=1:nS
            dumD(s,1:kd_s(s)) = 1;
        end
        
        % iterate on debt
        for iter_D=1:5
            for s=1:nS
                fun_lamD       = griddedInterpolant(k,lamD(s,:),meth,ex_meth);
                lamDp          = fun_lamD(kp_Q);
                Z_D(:,s)       = (exp_gp_Q.*lamDp)*probE;
                lamDp          = fun_lamD(kp_Q_iss);
                Z_D_iss(:,s)   = (exp_gp_Q_iss.*lamDp)*probE;
            end
            EMD     = reshape( sum(Z_D     .* Q_big,2) , [nS,nK] )./Rf1_mat;
            EMD_iss =          sum(Z_D_iss .* Q    ,2)            ./Rf(:,1); % [nS,1]
            lamDnew =           dumD.*(lamA.*(1-omega)) + (1-dumD-dumI).*((1-tau_i)./k' + EMD) + dumI .*((1-tau_i)./k' + (k_bar./k').*EMD_iss);
            dD      = max(abs(lamDnew(:)-lamD(:)));
            lamD    = lamDnew;
        end
        lamDx = lamD - bsxfun(@rdivide,1-dumD,k');
        
    end
    
    % optimal interest coverage ratio = k_bar
    lamF = lamS + lamD;
    for s=1:nS
        fun_lamF = griddedInterpolant(k,lamF(s,:),meth,ex_meth); % here 's' is tomorrow's state
        lamFp    = fun_lamF(kp_Q);
        Z_F(:,s) = (exp_gp_Q.*lamFp)*probE;
    end
    EMF         = reshape( sum(Z_F .* Q_big,2) , [nS,nK] )./Rf1_mat; % [nS,nK]
    noISS       = (1-tau_c)*(1-tau_d)*(1-1./k').*(1-equity_iss) + (1-tau_i)./k' + (1-tau_c+xi_e)*(1-1./k').*equity_iss + EMF; % [nS,nK] no issuance
    max_this    = EMF - xi_d*lamDx;
    [opt,ix]    = max(max_this,[],2); % find optimal kappa (on grid)
    ISS         = (1-tau_c)*(1-tau_d)*(1-1./k') + (1-tau_i)./k' + opt; % issuance (assumes: debt and equity issue do not occur simultaneously)
    k_bar_new   = k(ix); % [nS,1]
    dk_bar      = max(abs(k_bar_new(:)-k_bar(:)));
    k_bar       = k_bar_new;
    
    % optimal issuance threshold = k_iss
    dumI        = (ISS>noISS) & (k'>k_bar);
    ki_s        = min(nK,sum(1-dumI,2)+1); % index of k-value that equals issuance threshold
    k_iss_new   = k(ki_s); % [nS,1]
    
    dk_iss      = max(abs(k_iss_new(:)-k_iss(:)));
    k_iss       = k_iss_new;
    dumI        = k'>k_iss;
    
    %fprintf('Iteration: %4.0f, dk_bar = %4.6f, dk_iss = %4.6f, dS = %4.6f, dD = %4.6f \n',iter,dk_bar,dk_iss,dS,dD)
    %fprintf('Iteration: %4.0f, k_bar = %3.4f,%3.4f,%3.4f,  k_iss = %3.4f,%3.4f,%3.4f,  dS = %4.6f, dD = %4.6f \n',iter,k_bar(:),k_iss(:),dS,dD)
    
    if dk_bar==0 && dk_iss==0 && dS<prec && dD<prec && iter>20, break; end
        
end

kd_s    = sum(lamS==0,2);
k_def   = k(kd_s); % [nS,1]

% ex-dividend S/E ratio
lamSx    = (1-dumI).*EMS + dumI.*EMS_iss;
lamSx(lamS==0) = 0;
Div     = lamS - lamSx; % total dividend
SDiv    = dumI.*(1-xi_d-k_bar./k').*lam_Dx_kbar; % special dividend
RDiv    = Div - SDiv; % regular dividend

lamDx   = (1-dumI).*EMD + dumI.*(k_bar./k').*EMD_iss;
Coup    = lamD - lamDx;% total coupon


